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Overview  and  Research  Summary 


Introduction 

This  research,  funded  by  a  grant  from  the  Army  Research  Office  of  Terrestrial  Science 
Hydrology/Geomorpholgy  Program,  is  developing  a  physically-based,  Low-Dimensional  Modeling 
System  (LDMS)  for  runoff  for  regions  of  complex  terrain.  A  traditional  approach  to  rainfall-runoff 
focuses  on  channel  and  surface  flow  processes.  However,  this  study  emphasizes  the  longer  time 
scales  of  soil  moisture  and  subsurface  storage  in  controlling  runoff,  within  regions  of  complex 
topography  and  geology.  For  a  single  stream-reach,  our  modeling  strategy  identifies  how  hillslope 
soil  moisture  and  subsurface  storage  interact  with  channel  processes  through  direct  coupling  and 
feedback  from  baseflow  and  bank  storage,  deep  recharge,  and  saturation-overland  flow  processes. 
The  model  includes  the  essential  components  of  atmospheric  forcing  from  rainfall,  evaporation, 
transpiration,  snow  accumulation  and  snowmelt  and  their  dynamic  interaction  with  soil  moisture  and 
saturated  groundwater  to  produce  runoff.  In  the  case  of  stream  networks  and/or  vertically  stratified 
geology  contributing  to  runoff,  the  hillslope-stream-reach  becomes  an  element  in  an  extended 
dynamical  system.  The  extended  dynamical  system  is  a  spatially  distributed  model  which  attempts  to 
preserve  the  statistically  important  time/space  scales  of  the  catchment  runoff.  In  this  sense,  a 
particular  hillslope  and  corresponding  stream-reach  represents  a  spatial  random  variate  (Duffy, 

1996)  from  the  underlying  distribution  of  hillslopes  within  the  watershed  or  drainage  basin.  At  any 
point  within  the  stream  network,  runoff  represents  an  integration  of  the  upland  contribution  from 
spatially  random  elements  of  hillslope  storage  and  flux.  From  this  perspective,  the  question  for 
prediction  is:  What  is  the  distribution  of  contributing  hillslopes  and  what  are  their  "expected" 
time/space  scales?  From  our  signal  processing  and  dimension  estimation  work,  the  number  of  time 
scales  contributing  to  runoff  seems  to  be  limited,  with  only  a  few  hillslope  scales  contributing  most 
of  the  runoff.  The  number  of  states  contributing  to  runoff  becomes  a  problem  of  state-space 
dimension  estimation. 

This  report  will  outline  a  systematic  approach  to  the  problem  of  physically-based  modeling 
in  complex  terrain  involving  four  basic  steps: 

•  Dynamical  model  development  at  the  hillslope-stream-reach  scale. 

•  State-space  dimension  estimation  from  historical  streamflow  and  hydroclimatic 
observations. 

•  Conceptual-mathematical  model  development  based  on  terrain  features,  geology  and  the 
relation  to  the  stream  network. 

•  Model  and  parameter  identification/calibration  and  verification 
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The  first  phase  of  this  research  effort  has  involved  the  development  of  a  hillslope  dynamical 
modeling  system  which  couples  the  dynamics  (feedback  and  interaction)  of  the  important  physical 
processes  within  a  stream  reach.  A  paper  by  Duffy  (1996)  and  a  dissertation  by  Lee  (1994)  provides 
details  of  the  first  development  of  the  model.  Recent  progress  has  involved  examining  the  stability  of 
the  hillslope-stream-reach  system,  and  an  extension  of  the  model  to  include  the  effects  of 
macroporous  soils  and  rocks  on  soil  moisture,  groundwater  storage  and  runoff  timing.  A  paper  on 
the  stability  of  the  hillslope  component  of  the  model  was  recently  published  (Brandes,  Duffy  and 
Cusumano,  1998),  and  another  paper  has  been  completed  on  the  macropore  extension  of  the  model 
Brandes  and  Duffy,  2000). 

The  problem  of  dimension  estimation  and  the  forming  of  the  hydrogeologic  conceptual 
model  is  an  essential  ingredient  in  the  LDMS  approach.  Qualitatively,  the  characteristic  terrain  of 
physiographic  regions  of  North  America  implies  some  measure  of  similarity  in  form.  Hydrologists 
have  also  known  that  there  seem  to  be  a  limited  set  of  time/space  scales  contributing  to  runoff 
(Jakeman  and  Homberger,  1984).  Here  we  attempt  to  move  beyond  this  qualitative  notion  of 
"characteristic  terrain",  to  a  quantitative  statistical  method  to  estimate  the  dimension  of  the  state- 
space  contributing  to  runoff  at  a  point.  It  can  be  thought  of  as  the  number  of  independent  hillslopes 
(time/space  scales),  which  make  up  the  dominant  contribution  to  the  runoff  time  series.  We  have 
adopted  the  signal  processing  technique  of  singular  spectrum  analysis,  also  known  as  Karhunen- 
Loeve  decomposition  or  principal  orthogonal  decomposition  as  a  means  of  estimating  this  dimension 
from  historical  precipitation-temperature  runoff  data.  We  have  applied  the  approach  to  mountain- 
front  runoff  response  for  the  Wasatch  Front  (Shun  and  Duffy,  1998),  the  Susquehanna,  Rio  Grande, 
Colorado  and  Washita  rivers  in  the  US.  The  dimension  estimation  amounts  to  a  partitioning  of  time 
variance  of  runoff  into  independent  components.  A  ranking  of  the  eigenvalues  for  each  independent 
component  gives  a  measure  of  the  relative  contribution  of  each  independent  mode  (eigenvector)  to 
the  overall  runoff  signal.  The  dominant  modes  of  fluctuation  are  then  qualitatively  related  to 
characteristic  terrain  features  ,  soils  and  geology.  This  leads  to  a  conceptual  model  of  the  hydrology- 
hydrogeology,  the  hillslope-stream-reach  dynamcial  model  and  finally  to  an  estimate  of  the  extended 
dynamical  system  for  watershed. 

The  final  phase  of  the  research  is  the  development  of  a  calibration  and  validation  procedure  with 
field  applications  over  a  range  of  physiographic,  climatic  and  geologic  conditions  across  North 
America.  We  have  successfully  incorporated  an  optimization  algorithm  based  on  the  genetic 
algorithm.  The  genetic  algorithm  serves  two  purposes: 

•  The  algorithm  performs  a  random  search  of  the  expected  range  of  the  physical  parameter 
space  and  allows  the  user  to  utilize  available  hydraulic  and  physiographic  information  on  the 
soil  and  geologic  material  (e.g.  porosity,  hydraulic  conductivity),  terrain  features  (slope, 
length  of  flow,  contributing  area)  from  the  watershed  database. 
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•  It  allows  the  modeler  to  test  the  relative  significance  of  adding  new  state  variables  (hillslope- 
reach,  or  deeper  geologic  strata)  to  improve  the  model  performance. 

The  case  study  presented  here  is  the  Upper  West  Branch  of  the  Susquehanna  River.  This  paper  has 
recently  been  submitted  to  Water  Resources  Research  (Shun  and  Duffy,  2000).  A  second  model 
application  is  underway  in  the  Washita  River  basin  in  Oklahoma  and  this  paper  is  in  progress.  The 
third  application  is  the  Shale  Hills  watershed  near  State  College,  PA.  The  model  has  been 
implemented,  parameters  estimated  and  a  paper  will  be  completed  during  2000. 


The  Dynamical  Hillslope-Stream-Reach  Model 


A  dynamical  model  is  devised  for  a  where  unsaturated  and  saturated  storage  serve  as  the  principal 
control  on  rainfall-runoff,  and  where  complex  topography,  drainage  area  and  variable  depth  of 
moisture  penetration  describe  the  flow  geometry.  The  model  is  formed  by  direct  integration  of  the 
local  conservation  equation  with  respect  to  the  partial  volumes  occupied  by  unsaturated  and 
saturated  moisture  storage,  respectively  (Duffy,  1996).  This  yields  an  "integral-balance"  m.  The 
parametric  form  of  the  storage-flux  or  constitutive  relationships  for  the  proposed  model  is 
determined  from  numerical  experiments  in  a  simple  hillslope  flow  geometry. 

It  can  be  said  that  the  fundamental  equation  of  hydrology  is  the  water  balance  or  conservation  of 
fluid  volume 


dV 

dr 


-l-Q 


(1) 


where  V  represents  the  global  fluid  storage  volume,  1  the  volumetric  input  and  Q  the  volumetric 
discharge  or  output.  For  a  hydrologic  system  with  interacting  physical  processes,  (1)  can  be  written 
in  state-space  form: 


dS 

dt 


=f-  g 


(2) 


where  the  state  variables  or  storage  components  form  a  vector  S  ={Sl,S2,...S^};  the 
vector  s/  =  {/,,/2, and  g=  {gi,g2f-g„}  represent  the  various  input-output  flux  components 
to  and  from  each  state.  Both  fandg  are  in  general,  functions  of  the  system  states,  5. 

Although  the  state-space  framework  offers  a  simple  representation  of  the  system  dynamics,  the 
particular  form  and  meaning  of  the  inputs,  outputs  and  states  must  be  made  explicit.  The 
interpretation  of  (2)  as  a  conservation  statement  for  the  system  implies  that  the  state  variables  are 
formed  by  integration  over  an  appropriate  flow  volume  for  each  state,  and  that  the  observables  are 
defined  as  volume  averages.  Likewise,  the  flux  quantities  represent  integrations  of  the  distributed 
flow  over  the  surfaces  separating  the  states.  The  relationship  between  the  surface  fluxes  and  states 
represent  the  constitutive  relations  for  the  system.  In  the  next  section  an  integral-balance 
formulation  is  developed  at  the  hillslope  scale  which  represents  the  elemental  spatial  scale  for 
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predicting  runoff  in  complex  terrain.  Besides  developing  the  basic  model  equations  this  section 
outlines  the  physical  processes  which  relate  storage  to  flux  and  define  the  constitutive  relations  of 
the  model. 


The  integral-balance  for  the  hillslope-stream-reach  in  state-space  form  can  be  written: 


d5j 
d  t 


f  01  Sl2  J  1 


10 


dt 


—  f  02  &  1  &3  f 


20 


d  & 

TT  =  <?03  +  #23  “  <ho 


(3) 


where 


S',  =  unsaturated  moisture  storage 

52  =  saturated  moisture  storage 

53  =  channel  storage 

/0,  =  incoming  boundary  flux  to  5, 

/l0  =  outgoing  boundary  flux  from  S, 

f20  =  outgoing  boundary  flux  from  S2 

f  Q2  =  incoming  boundary  flux' to  S2 

g2 3  =  -g31  =  groundwater  -  channel  flux  S2  <=>  S^ 

g21  =  -gl2  =  recharge  flux  S <=>  S2 

qa 3  =  upstream  inflow  to  channel  reach 

q30  =  downstream  outflow  from  reach 

The  double  subscripts  on  boundary  flux  quantities  define  the  predominant  sense  of  the  flow  direction 
as  from-to  (O-exteranl  environment,  1-unsaturated  storage,  2-saturated  storage,  3-channel  storage). 
For  example  the  subscript  01  refers  to  a  flow  from  the  environment  directed  to  state  1  (infiltration). 
The  next  section  deals  with  the  relation  of  state  variables  to  flux  rates  for  the  hillslope  system  and 
the  numerical  experiments  carried  out  to  define  each  relationship. 

Numerical  Experiments  and  Development  of  Flux-Storage  Relationships 

The  integration  of  the  local  equations  removes  the  infinite  dimensionality  associated  with 
continuously  varying  parameters  and  variables  in  space  (as  in  the  PDEs),  and  also  provides  a 
hydrologically  useful  scale  for  modeling  -  that  of  the  hillslope.  The  effects  of  local  processes  and 
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local  spatial  heterogeneity  on  integral  hillslope  response  must  be  parameterized  through  the  integral- 
scale  constitutive  relations,  that  is,  “effective”  representations  of  such  processes  are  required.  The 
ability  to  simulate  at  spatial  scales  less  than  that  of  the  hillslope  is  obviously  lost;  however,  as  a 
practical  matter  such  detail  is  usually  unnecessary  in  modeling  fluxes  relevant  to  water  resource 
management. 

To  develop  the  model,  one  needs  to  determine  how  the  integrated  moisture  state  variables  are 
related,  i.e.,  the  equilibria,  for  typical  hillslope  geometries.  In  addition,  the  form  of  the  integral 
constitutive  relations  for  hillslope  fluxes  f(SvS2)  and  g(Sp5,)  must  be  developed.  Given  such 
relations,  the  system  (3)  becomes  a  low-dimensional,  physically-based  model  of  hillslope-stream- 
reach  hydrology. 

The  key  to  formulating  a  useful  model  from  the  integral-balance  structure,  is  the  determination 
of  appropriate  hillslope  or  catchment-scale  constitutive  relations.  By  incorporating  these  relations 
into  the  integral-balance,  the  model  becomes  system-specific.  The  approach  we  have  taken 
(Lee, 1993;  Duffy,  1996;  Brandes,  1998)  is  to  base  these  relations  upon  steady-state  numerical 
solutions  of  Richards’  equation  over  a  hillslope  domain.  Implicit  in  this  approach  is  the  assumption 
that  the  dynamics  far  from  equilibrium  are  in  some  sense  similar  to  those  near  equilibrium,  and  can 
be  represented  by  the  steady-state  storage-flux  relations 

Approach 

Richards’  equation  in  pressure-head  form  (1)  is  solved  over  a  two-dimensional  discretized  hillslope 
domain  using  the  galerkin  finite  element  code  FEMWATER  (Yeh,  1987).  The  nonlinear  matrix  of 
element  equations  is  linearized  by  guessing  initial  values  for  {y/},  and  then  is  solved  directly  by 

Gaussian  elimination.  New  values  of  {i //}  are  estimated  using  under- relaxation,  and  the  procedure  is 
repeated  until  convergence  (i//  to  within  0.0001  meter)  is  achieved.  For  approximating  the  time 
derivative,  the  fully  implicit  method  was  used.  Further  details  of  the  numerical  method  are 
contained  in  Lee  (1993). 

The  experiments  involve  steady-state  solutions  for  a  range  of  constant  precipitation  rates  to  cover  the 
full  range  of  hillslope  saturation.  Precipitation  is  applied  uniformly  to  the  hillslope.  The  numerical 
runs  were  made  from  an  arbitrary  initial  condition,  usually  a  constant  pressure  head  throughout  the 
domain.  For  each  experiment,  the  integrated  moisture  state  variables  S,  and  S2  are  determined 
numerically  by  summing  over  the  saturated  and  unsaturated  hillslope  volumes.  The  steady-state 
position  of  the  water  table  and  associated  surface  saturated  area  are  also  determined  during  the 
experiment.  In  addition,  the  code  provides  total  flux  rates  through  the  various  boundaries  (stream, 
seepage,  infiltration).  Note  that  at  steady-state,  the  subsurface  flux  rate  will  generally  not  equal  the 
precipitation  rate  applied,  due  to  rejection  of  precipitation  along  the  saturated  portion  of  the  ground 
surface.  This  physical  process,  known  as  saturation  overland  flow,  is  an  important  outcome  of  the 
steady-state  numerical  experiments. 
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Hillslope  geometry,  boundary  conditions  and  soil  properties 
The  geometry  used  for  most  of  our  experiments  is  that  of  a  convex-concave  hillslope  of  height  (H) 
25  m,  length  ( L )  100  m,  and  soil  depth  ( d )  of  2.5  m  over  an  impermeable  base  (see  Figure  1).  These 
shape  function  is  considered  typical  of  hillslopes  in  humid-temperate  regions.  The  surface  shape 
function  is  given  by  (Lee,  1993): 


h(x)  = 


H 


e4  -1 


(4) 


rTTTTTTTTTTTTTl 


P 


Figure  1.  Schematic  of  the  hillslope  geometry  and  boundary  conditions  used  in  the  numerical 
experiments  (not  to  scale). 


The  domain  is  discretized  into  4000  quadrilateral  elements  (4221  nodes),  consisting  of  200  rows  by 
20  columns.  Thus,  each  element  is  0.5  m  long  by  0.125  m  high. 

The  boundary  conditions  include  no-flow  (Neumann)  boundaries  along  the  sides  and  base  of  the 
hillslope,  a  variable  infiltration/seepage  boundary  along  the  ground  surface,  and  a  single  constant 
head  (Dirichlet)  node  at  the  foot  of  the  slope  representing  a  first-order  stream.  The  no-flow  boundary 
beneath  the  stream  is  due  to  converging  flow  from  the  hillslopes  on  opposite  sides  of  the  stream 
(symmetry  condition),  and  the  no-flow  boundary  along  the  base  of  the  hillslope  represents  an 
impermeable  bedrock  layer. 
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For  this  report,  results  for  three  soil  types  are  reported,  a  Guelph  loam,  a  Plainfield  sand,  and  a 
medium  sand  (Clapp  and  Homberger,  1977)  designated  Cl&H  sand.  These  soils  represent  a  wide 
range  of  hydraulic  characteristics.  Additional  soil  types  and  hillslope  geometries  have  been 
investigated  by  Lee  (1993)  and  Brandes  (1998),  and  those  presented  here  reflect  the  general  results. 
The  numerical  code  requires  input  of  analytical  functions  or  tabular  data  describing  soil  moisture 
content,  soil  moisture  capacity,  and  unsaturated  hydraulic  conductivity  as  a  function  of  pore 
pressure.  Analytical  soil  hydraulic  property  functions  were  generally  used  in  the  numerical  model; 
however,  in  some  cases  (very  low  precipitation  rates,  highly  permeable  soils)  tabular  input  was 
necessary  for  numerical  convergence.  For  the  Guelph  loam  and  Plainfield  sand,  soil  hydraulic 
properties  were  described  by  the  nonhysteretic  Elrick  et  al.  (1990)  parameterizations: 

W=  +  c<  (5) 

c2+m 

K(y/)  =  Ksate-k*  (6) 

where  6 ,  ly ,  and  K  are  defined  as  previously,  Ksal  is  the  saturated  hydraulic  conductivity,  and  ci  and 

k  are  Fitting  parameters.  For  consistency,  the  Cl&H  sand  property  data  were  also  fitted  to  these 
parameterizations  rather  than  the  method  used  by  Clapp  and  Homberger  (1977).  Table  1 
summarizes  the  soil  property  parameter  values  used  in  the  numerical  experiments: 


Soil  Type 

Ksal 

(cm/s) 

es 

C, 

c2 

c3 

C4 

Guelph  loam 

3.67e-04 

811111 

3.36 

0.243 

0.421 

2.0 

0.28 

Plainfield 

sand 

3.44e-03 

13.06 

0.377 

0.00154 

0.1 

Cl&H  sand 

1.76e-02 

0.395 

6.24 

0.323 

0.927 

0.086 

Table  1  Summary  of  soil  property  data  used  in  the  numerical  experiments. 

The  selected  soils  cover  a  range  of  almost  three  orders  of  magnitude  in  conductivity  and  a  wide 
range  in  porosity.  The  water  retention  and  relative  conductivity  functions  are  shown  in  Figure  2. 
Note  that  the  Guelph  loam  retains  significantly  more  moisture  than  the  sands,  and  thus  shows  a 
gradual  decrease  in  conductivity  with  pore  tension,  while  the  Plainfield  sand  exhibits  the 
characteristics  of  a  uniformly-graded  soil,  with  moisture  content  and  conductivity  dropping  sharply 
over  a  narrow  range  of  pore  tension.  The  Cl&H  sand  hydraulic  behavior  is  intermediate  to  the 
Guelph  and  Plainfield  soils.  Note  that  channel  stage  effects  are  included  in  these  experiments  and 
thus  the  emphasis  is  on  developing  relationships  for  the  soil  moisture  and  saturated  storage. 
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Figure  2  a)  Soil  moisture  characteristic  curves,  and  b)  hydraulic  conductivity  functions  for  the 
Guelph  loam,  Plainfield  sand,  and  Cl&H  sand  ( Kret  =  K(y)/Ksat). 


Results 

Results  for  the  Guelph  loam  soil  are  shown  in  Figure  3.  The  first  plot  ( a)  shows  the  steady-state 
relation  St*-S2*  ,  developed  over  a  range  of  precipitation  rates,  i.e.,  the  fixed  points  for  different 
values  of  precipitation  p.  Note  that  the  *  indicates  that  the  depth  of  storage  is  made  dimensionless  by 
the 'd'  the  depth  of  the  soil.  The  *  is  dropped  in  later  plots.  The  second  plot  (b)  shows  S,  and  S2 
against  f0I,  the  net  steady  flux  through  the  hillslope.  The  third  plot  (part  c)  shows  fractional  saturated 
surface  area  (frSSA )  against  S2.  Similar  results  were  obtained  for  the  other  soil  types. 

A  fundamental  result  of  these  numerical  experiments  is  the  competitive  or  inverse 
relationship  between  the  moisture  state  variables  over  most  of  the  range  of  soil  moisture,  except 
when  the  water  table  is  at  or  near  its  lowest  state  (see  Figure  3a).  The  form  of  this  relationship  can 
be  understood  in  the  following  way:  the  hillslope  contains  a  fixed  volume  of  pore  space,  that  must 
be  shared  by  the  saturated  moisture  ( S2 ),  the  unsaturated  moisture  ( S, ),  and  the  empty  (air-filled)  pore 
space.  Thus,  except  where  the  empty  pore  volume  is  large,  there  is  an  inverse  or  competitive 
relationship  between  S,  and  S2.  This  general  structure  is  seen  for  a  variety  of  soil  types  and  hillslope 
configurations  (Lee,  1993;  Brandes,  unpublished).  Based  on  extensive  numerical  experiments,  the 
flux-storage  or  integral  constitutive  relations  are  found  to  be  nonlinear. 
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Figure  3  Numerical  experiment  results  for  the  Guelph  loam  hillslope:  a)  fixed  points  fir  ,  b)  flux- 
storage  relations,  and  c)  saturated  surface  area-storage  relation.  f0l  is  the  steady  flux  rate  through  the 
hillslope;  frSSA  is  the  fractional  surface  saturated  area.  Fits  to  the  model  are  shown  as  solid  lines. 


The  fitted  two-state  variable  model 

Duffy  (1996)  parameterized  the  results  shown  in  Figure  3  above,  and  developed  the  following  forms 
for  the  flux  terms  of  for  a  two-state  model  (no  channel  storage): 


foi 

=  p[l-djs2-s°2)] 

(7) 

Sl2 

=  d^S^Hd^+d^S.-Sl)2 

(8) 

ho 

=  d3(S2-S°2) 

(9) 

in 


The  form  of  the  model  becomes: 

^■  =  p[>-  -  d,(S,  +  d2)(S2-S°2)\  10a) 

^  =  d0(s,-s;>+  d,(S,  +d2)(S2~S2)1  -  d1(S,-S‘2)(  10b) 

at 

where  p  is  the  precipitation  rate,  S°  and  S%  represent  residual  storage  {p  =  0),  and  the  d’s  are  fitted 
parameters.  The  term  f0l  is  the  net  infiltration  to  the  hillslope.  This  “effective  precipitation”  is 
inversely  proportional  to  the  extent  of  surface  saturation,  which  is  itself  a  linear  function  of  saturated 
storage.  The  rejected  precipitation  p^ip  -  f0I)  becomes  ssurface  runoff  to  the  channel,  and  thus  the 
model  represents  the  mechanism  of  saturation-excess  runoff  in  addition  to  subsurface  flow.  Outflow 
from  the  hillslope  (f20 )  is  shown  to  be  a  function  of  saturated  storage  (linear  in  eq.9).  The  steady- 
state  rescharge  rate  gI2  contains  a  nonlinear  coupling  for  the  two  state  variables.  The  justification  for 
this  form  of  the  model  is  that  it  is  consistent  with  the  numerical  experiment  results  and  is  of  low 
order  (Duffy,  1996);  however,  other  forms  of  the  model  are  certainly  possible,  for  instance  quadratic 
or  cubic  external  flux  relations.  The  fluxes  fI0  (evaporation)  and  fm  (leakage  into  the  saturated  zone) 
were  not  included  in  these  numerical  results  and  will  be  developed  in  a  later  section. 

The  physical  meaning  of  the  model  parameters  can  be  generally  deduced  from  the  flux  terms  (7)  - 
(9).  These  parameters  can  be  thought  of  as  integrated  versions  of  the  parameters  describing  the  C( if/) 

and  K( ip)  relationships  of  Richards’  equation  at  the  local  scale.  The  parameters  d0,  d„  and  d2 

describe  the  form  of  the  moisture  relation  5;(52),  which  is  a  function  of  soil  hydraulic  properties  and 
hillslope  geometry.  Parameter  d}  is  a  rate  constant  for  subsurface  flow,  which  depends  on  hillslope 
geometry  and  soil  hydraulic  conductivity.  Parameter  d4  describes  the  relation  between  surface 
saturation  and  saturated  storage,  a  function  of  hillslope  geometry. 

In  terms  of  the  saturated  zone,  this  physically-based  formulation  is  not  unlike  TOPMODEL  of 
Beven  and  Kirkby  (1979),  which  includes  the  processes  of  infiltration  and  saturated  storage,  surface 
runoff  from  saturated  areas  (quick  flow  response),  and  subsurface  runoff  (delayed  response). 
However,  in  their  model,  infiltration  simply  cascades  to  the  saturated  zone  when  a  threshold  is 
reached.  In  contrast  with  (7)  the  competitive  and  nonlinear  cross  term  for  recharge  is  a  critical 
component  of  the  model,  providing  two-way  coupling  between  the  unsaturated  and  saturated 
dynamics  consistent  with  the  physics  of  recharge  and  relaxation  of  the  water  table  surface  (Duffy, 
1996). 

Macropore  Flow  and  Hillslope  Dynamics 

Numerical  experiments  carried  out  by  Brandes  (1998)  and  Lee  (1993)  for  unsaturated- 
saturated  flow  on  hillslopes  using  classical  parameterizations  of  matrix  soil  hydraulic  properties, 
indicate  excessively  long  response/relaxation  times  in  comparison  to  similar  field  soils.  Apparently 
the  hydraulic  properties  of  field  soils  are  much  smaller  than  the  integrated  effective  conductivity  of 
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the  hillslope  or  watershed.  We  attribute  this  significant  "error"  to  "macropore"  flow  in  the  field  soils, 
that  is,  flow  through  interpedal  cracks,  root  holes,  or  other  channels  of  biological  origin  not 
accounted  for  in  the  soil  matrix  hydraulic  conductivity.  Further,  recent  data  for  macroporous  soils 
suggest  hydraulic  conductivity  functions  ( K( if/))  often  have  a  sharp  increase  in  conductivity  near 

saturation  and  an  overall  bimodal  appearance  of  the  K(yf)  relationship. 

A  simple  method  is  devised  for  estimating  such  hydraulic  conductivity  functions  for 
macroporous  soils,  conceptualizing  the  soil  as  an  equivalent  porous  medium  with  pore-size-weighted 
effective  conductivity.  The  method  requires  a  macropore  size  density  function  and  macroporosity 
estimate,  in  addition  to  the  matrix  hydraulic  conductivity  function  and  total  porosity.  The  resulting 
effective  K( yr)  functions  are  then  used  to  numerically  model  variably  saturated  flow  on  a  hillslope. 

The  numerical  experiments  suggest  that  the  presence  of  macropores:  1)  increases  the  steady-state 
moisture  storage  in  the  unsaturated  state;  2)  increases  the  steady  outflow  rates  from  saturated  storage 
in  proportion  to  the  increase  in  effective  saturated  conductivity;  3)  decreases  the  subsurface  storm 
response  time  (time  to  peak);  and  4)  significantly  alters  the  shape  of  the  storm  hydrograph  and  the 
relative  contributions  of  surface  and  subsurface  flow.  A  brief  discussion  of  the  method  and  results  is 
presented  next. 

“Macropore”  flow  through  litter,  structured  soils  or  biologically-induced  channels  is  quite 
common  on  forested  slopes  (Beven  and  Germann,  1982).  Observations  of  rapid  subsurface  response 
at  many  humid  forested  watersheds  (see  for  example  Hursh,  1936;  Hursh,  1944;  Whipkey,  1965)  and 
a  semiarid  hillslope  site  (Wilcox  et  al.,  1997)  suggest  that  macropore  flow  is  an  essential  component 
of  the  runoff  response  of  upland  watersheds.  Moreover,  fitting  of  recession  flow  data  from  small 
watersheds  to  lumped-parameter  hydraulic  models  (by  inverse  methods)  has  resulted  in  basin-scale 
effective  saturated  conductivity  values  approximately  one  to  two  orders  of  magnitude  greater  than 
indicated  by  available  laboratory  measurements  (Zecharias  and  Brutsaert,  1988;  Troch  et.  al.,  1993). 

We  first  develop  a  simple,  yet  physically-based  parameterization  of  the  hydraulic 
conductivity  function  for  macroporous  soils.  The  approach  taken  here  is  to  modify  the  matrix  K( y/) 

function  to  include  the  effect  of  macropore  flow  near  saturation.  Since  the  focus  is  on  hillslope-scale 
behavior,  the  continuum  or  "equivalent  porous  medium"  approach  is  used  rather  than  treating  the 
macropore  domain  explicitly  (e.g.,  Barenblatt  et  al.,  1960;  Beven  and  Germann,  1981). 

Next  the  composite  or  effective  conductivity  functions  are  tested  using  numerical 
experiments  to  assess  the  effect  of  macropore  soil  properties  at  the  hillslope  storage  and  flow.  Of 
primary  interest  are  the  effects  on  moisture  equilibria,  outflow  time  scales,  surface  and  subsurface 
runoff  response.  The  numerical  experiments  include  steady-state  and  dynamic  conditions,  a  convex- 
concave  hillslope,  two  distinct  soil  types,  with  two  levels  of  macroporosity.  Results  are  compared  to 
those  for  the  same  soils  without  macropores. 

The  method  is  based  on  the  following  simplifying  assumptions: 
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•  During  conditions  when  macropore  flow  is  occurring  (soils  near  saturation),  the  pressure 
in  the  matrix  and  macropores  is  equal.  Therefore  transient  local  conditions,  such  as  water 
flowing  into  the  macropores  from  the  surface  and  then  diffusing  into  the  surrounding 
matrix,  are  not  explicitly  considered. 

•  The  above  assumption  implies  that  pressure  equilibration  occurs  rapidly  as  water  flows 
within  the  two  pore  domains. 

•  The  macropores  are  idealized  as  interconnected  capillary  tubes,  which  transmit  flow  only 
when  the  pore  tension  is  less  than  the  air  entry  suction  value.  Film  flow  along  macropore 
walls  and  the  potential  effects  of  pore  channel  tortuosity  are  not  considered  here. 

Under  the  conditions  stated  above,  one  can  define  an  effective  hydraulic  conductivity  (Keff), 
which  is  a  porosity-weighted  average  of  the  matrix  conductivity  and  the  macropore  conductivity: 

Keff(Y)  -frhiP % mp(W )  +(1-  frMp)KMx(V)  (I*) 

where  fr^  is  the  fraction  of  the  total  porosity  composed  of  macropores,  KUP  is  the  conductivity  of 
the  macropore  network,  and  Km  is  the  conductivity  of  the  soil  matrix.  Both  KUP  and  A^will  be 
functions  of  the  pore  tension  yp  f°r  the  matrix  through  a  known  Kmiyf)  relationship  (as  determined 

in  the  laboratory),  and  for  the  macropores  through  capillarity,  which  will  govern  the  pore  sizes 
which  are  saturated  and  thus  transmitting  water  at  a  given  tension.  Principles  of  hydraulics  of  flow 
in  cylindrical  tubes  are  used  to  calculate  hydraulic  conductivity  of  the  macropores  (similar  to  the 
approach  used  by  Chen  and  Wagenet,  1992). 

To  apply  these  equations,  two  additional  pieces  of  data  are  required:  the  macroporosity 
fraction  frMP  and  a  macropore  pore  size  distribution  or  density  function/r).  Values  of  0, 0.001,  and 
0.01  were  used  for  frSiP  to  represent  a  full  range  of  reasonable  values.  Various  studies  have  found 
that  numbers  of  soil  pores  of  a  given  size  class  are  inversely  proportional  to  pore  size  (e.g.,  Wilson 
and  Luxmoore  (1988)  infiltrometer  data).  Therefore,  an  exponential  probability  density  function 
was  used: 


f(r)  =  k,exp(-k2r)  (12) 

where  r  is  the  pore  radius,  and  k,  and  k2  are  fitting  parameters.  Since  the  area  under  the  density 
function  f(r)  must  integrate  to  unity,  the  parameters  k,  and  k2  are  interdependent  (furthermore, 
approximately  equal  for  k ,  and  k2  >  1).  The  sensitivity  of  Keff  to  these  parameters  was  assessed,  and 
it  was  determined  that  the  effect  was  small  in  comparison  to  the  value  used  for  frMP.  For  the  results 
shown  here,  the  values  kt  =  0.63  and  k2  =  0.6  were  used  with  a  macropore  radius  range  of  0.01  to  2.5 


mm.  The  0.01  mm  radius  corresponds  to  a  if/^  value  of  -1.5  m;  below  this  value,  flow  is  assumed  to 
occur  only  in  the  matrix  domain. 

Given  the  function/;*),  the  weighted  macropore  conductivity  is  calculated  using  the 
following  equation: 

r* 

j  KMP(r  )f(r)dr 

KmpW  -  <13> 

/  f(r)dr 
0 

with  r*  chosen  such  that  if/ae(r*)=\f/,  that  is,  the  weighted  sum  includes  only  those  pores  (with  r  <  r*) 
which  are  contributing  to  flow  at  a  given  value  of  iff.  This  integral  is  evaluated  discretely  for 
increasing  pore  tensions  (decreasing  values  of  r*),  starting  from  saturated  conditions  where  all  pores 
are  contributing  flow. 


Matrix  soil  properties 

Two  soil  types  were  used  in  the  Keff  calculations  and  simulations,  a  Guelph  loam  and  a 
Plainfield  sand.  The  Guelph  loam  has  finer-grain  size,  and  retains  significantly  more  moisture  than 
the  sand,  and  thus  shows  a  gradual  decrease  in  conductivity  with  pore  tension.  The  Plainfield  sand 
exhibits  the  characteristics  of  a  uniformly-graded  soil,  with  moisture  content  and  conductivity 
dropping  sharply  over  a  narrow  range  of  pore  tension  between  -0. 1  m  and  -0.5  m.  The  Km( iff) 
functions  for  the  unsaturated  soils  is  represented  by  the  exponential  Gardner  relation: 

KMX (¥)=  Ksa:  exp(-a Iff)  ( 14) 

where  Ksat  =  3.67X10-4  cm/sec  and  a  =  3.36  m'1  for  the  Guelph  loam,  and  Ksal  -  3.44xl0'3  cm/sec  and 
a=  13.06  m'1  for  the  Plainfield  sand  (Elrick,  1990).  These  two  soils  bracket  the  threshold  saturated 

conductivity  value  of  10‘3  cm/sec  for  subsurface  storm  flow,  suggested  by  Freeze  (1972). 

Once  the  composite  conductivity  functions  were  developed,  the  galerkin  finite  element  code 
FEMWATER  [Yeh,  1987]  was  again  used  to  solve  Richard's  equation  over  a  two-dimensional 
hillslope  domain  as  described  earlier. 

Steady-state  results  for  the  two  soil  types  and  three  values  of  macropore  fractions  ( f  - 0.0, 
0.01  and  0.001)  are  shown  in  Figures  4  through  7.  Figures  4  and  5  show  the  steady-state  integrated 
moisture  ( S2 *,  S,*)  and  storage-flux  ( q(S2 ))  for  a  range  of  constant  precipitation  rates,  and  shows  the 
expected  inverse  relationship  between  saturated  and  unsaturated  moisture,  except  under  dry 
conditions,  as  discussed  previously  by  Duffy  (1996)  and  Brandes  et  al.  (1998).  We  note  that 
although  the  soil  moisture  retention  curve  6(if/)  is  unchanged  by  the  macropore  parameterization, 

there  is  a  significant  effect  on  the  unsaturated-saturated  storage,  since  this  depends  on  the  overall 

14 


rate  of  flux  through  the  hillslope.  In  particular,  for  macropore  flow,  higher  flux  rates  and  a  lower 
water  table  are  observed  within  the  saturated  zone  to  balance  the  increased  steady  flux  rate  through 
the  hillslope  .  The  lower  water  table  leads  to  an  overall  higher  moisture  storage  in  the  unsaturated 
zone. 

The  effect  of  macropore  flow  is  most  pronounced  for  the  Plainfield  sand  (Figure  5).  For frMP-' 0.01, 
the  relation  attains  a  maximum  with  relatively  flat  slope  over  intermediate  values  of  S2 ■  In  this  case, 
the  macropore  conductivity  increases  over  the  same  range  of  pore  pressure  where  moisture  content 
changes  rapidly  in  the  soil  moisture  characteristic  6{y/).  Therefore,  relatively  large  increases  in  S, 
are  required  to  increase  the  unsaturated  hydraulic  conductivity  necessary  to  transmit  the  macropore- 
enhanced  fluxes. 

Under  very  dry  conditions  and  with  macropore  flow,  the  numerical  solutions  for  both  soils 
did  not  converge  near  the  peak  of  the  soil  moisture  storage  .  This  is  due  to  the  highly  nonlinear 
effective  conductivity  relations,  and  to  the  presence  of  weakly  damped  equilibria  (see  discussion  in 
Brandes  et  al.,  1998). 

The  steady  flux  data  are  shown  in  the  lower  portion  of  Figures  4  and  5,  give  the  subsurface 
flux  rate  as  a  function  of  the  saturated  storage  S2.  The  steady  flux  rates  for  the  Guelph  loam 
hillslope  increase  by  about  an  order  of  magnitude  for frMP=0.00l,  and  by  two  orders  of  magnitude  for 
frMP=0.0l.  The  effect  is  less  pronounced  for  the  Plainfield  sand.  Note  that  the  steady-state  fluxes  are 
not  equal  to  the  applied  precipitation  rate,  due  to  rejection  of  infiltration  along  the  saturated  surface 
near  the  stream.  This  relationship  between  the  extent  of  surface  saturation  and  saturated  storage  is 
determined  by  hillslope  geometry  and  boundary  conditions,  and  is  not  effected  by  macropores. 

It  is  clear  that  the  subsurface  flux  rates  are  greatly  enhanced  by  macropore  flow  as 
parameterized  here.  Interestingly,  this  effect  is  apparently  proportional  to  the  effective  saturated 
conductivity  of  the  hillslope.  Figure  6  shows  the  combined  data  from  Figures  6  and  7,  with  each 
curve  scaled  by  the  associated  Keff_3at  value.  Note  that  the  data  fall  very  close  to  a  single  curve, 
except  at  very  high  levels  of  saturated  storage  (Sp> 0.8).  Note  that  the  flux  -storage  relation  now 
requires  a  nonlinear  function  and  we  adopt  the  form: 

q  /  Ksat  =  [c,  (S2 -S°)+c2  (S2  -S°)2]  (15) 

where  c ,  =  0.0105,  c2  =  0.0032  ,  and  S2°  =  0.166  for  these  data.  The  c,  can  be  thought  of  as  rate 
constants,  and  the  parameter  S2°  describes  the  residual  moisture  storage  (corresponding  to  p=0). 
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Figure  4. .  Saturated-unsaturated  storage  (upper  plot)  and  the  subsurface  for  the  Plainfield  sand  for 
3  values  of  macropore  flow  frMP=0.0,frUP=0.0l  and Jr ^IP=0. 001,  (after  Brandes  and  Duffy,  2000). 
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Figure  5.  Saturated-unsaturated  storage  for  the  Guelph  loam  for  3  values  of  macropore  flow 
frMP=0.0,  frup=0.0 1  and  frMP=0.00U  (after  Brandes  and  Duffy,  2000). 
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Figure  6.  Scaled  data  for  both  soils  based  on  equa.  (14)  q  /  Ksat  =  f(S2)  (after  Brandes  and  Duffy, 
2000). 


Evaporation  and  Transpiration 

The  two-state  integral-balance  model  for  unsaturated-saturated  flow  (Duffy,  1996)  is 
extended  to  study  the  timing  of  atmosphere-controlled  and  soil-moisture  controlled  evaporation  for 
sloping  topography.  The  modeling  strategy  is  to  construct  the  state-space  for  unsaturated-saturated 
storage  by  direct  integration  of  the  unsaturated-saturated  flow  equations  (Richards'  equation)  and  to 
devise  a  parametric  model  for  evaporation.  An  empirical  evaporation  equation  is  proposed,  which  is 
a  function  of  both  the  unsaturated  and  saturated  state  variables.  The  dynamical  system  now  takes 
the  form 


=foi~8n  ~ew 

d S1_  _ 

it  “ U-  &l  30  /m  (  16 ) 

For  the  saturated  surface  (/3(S2))  evapotranspiration  flux  is  subject  to  "atmospheric"  control  and  e20 
is  simply  defined  as  the  product  of  the  potential  evaporation  pe and  ()3 (S2)): 

e20  =  PJ3(S2)  (17) 

Evapotranspiration  flux  over  the  unsaturated  surface  of  the  hillslope  is  assumed  to  be  controlled  by 
the  soil  moisture  and  is  defined  as  product  of  pe,  the  fraction  of  the  hillslope  that  is  unsaturated  and 
the  soil  evapotranspiration  function  w(^): 
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e10  =  ft(l-/3(S,))-H>tt) 


(18) 


The  soil  moisture  evapotranspiration  function  w(%)  is  a  weighting  function  (0,1)  defined  as: 


w(x)  = 


a+xm 


;  x  = 


s i 

d-n-  S2 


(19) 


The  function  x(Sl,S2)  represents  the  degree  of  saturation  of  the  unsaturated  volume,  d  the  average 
depth  of  the  soil  and  m  determines  the  shape  of  w(x)  •  The  formula  of  Thomthwaite  [1948]  relates 
potential  evapotranspiration  to  temperature  (ETp  °cTa)  where  T  is  the  mean  daily  air  temperature 
For  simplicity,  the  potential  ET  is  expressed  as  the  power  function  of  the  daily  maximum 
temperature, 


ETp  =  bx 


(  j  V2 

T 

V  1  max  J 


-b. 


(20) 


The  ET  model  results  suggest  that  the  extent  and  timing  of  surface-saturation  during  cycles  of 
wetting  and  drying  is  dependent  on  the  saturated  (storage)  of  the  watershed,  which  in  turn  is  closely 
tied  to  the  soil  moisture  storage  and  the  magnitude  of  surface  and  subsurface  runoff  components. 
The  transition  from  atmosphere-controlled  evaporation  to  soil-controlled  evaporation  following 
precipitation  is  examined  for  a  range  of  soil  types  and  hillslope  geometries.  This  transition  may  be 
abrupt  or  smooth  depending  on  the  parametric  form  of  the  evaporation  term  and  the  overall  storage 
characteristics  of  the  terrain. 


Snow  Accumulation  and  Melt 

A  simple  model  for  daily  potential  snowmelt  model  was  developed  by  the  U.S.  Army  Corps  of 
Engineers  and  described  by  Singh,  [1992]  and  have  adapted  the  model  here.  The  basic  idea  is  that 
the  potential  flux  due  to  snow  melting  is  a  simple  function  of  air  temperature 

M=0A{T-Tmelt)  (21) 

where  T  is  the  daily  maximum  air  temperature,  M  [mm/day]  is  the  potential  snowmelt  at  T  [°F]and 
T^,,  is  the  local  melt  temperature.  The  range  of  Tmelt  is  from  27  to  42  [°F]  which  corresponds  to 
open  and  forested  sites. 
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Based  on  the  potential  snowmelt  relation  (21),  the  accumulation  and  melt  of  snow  and  the  relation  to 
infiltration  is  formed  by  the  following  conditional  statements: 


Pm(t)  =  f{Ps,t)MT>Tmelt  (23) 

Pm(t)  =  0,  and  Ps (?)  =  Ps(t  - 1)  ■ +  P{t) ,  if  T  <  Tmell  (24) 

/(/>,?)  =  M  and  Ps(t  +  l)  =  Ps (?)  -  M,  if  M <  Ps (?)  (25) 

f(Ps>T)=  Ps(t)  if  M>Ps{t)  (26) 


where  Pm( t)  is  the  snowmelt  at  time  t,  P(t)  is  the  precipitation  at  time  t  and  Ps(t)  is  the  amount  of  the 
accumulated  snow-pack  at  time  t.  At  any  time,  precipitation  is  the  sum  of  rainfall  and  snowmelt. 


Space-Time  Decomposition  and  Dimension  Estimation 

Multichannel  SSA  is  a  space-time  version  of  the  Karhunen-Loeve  theorem  and  a 
general  discussion  of  the  method  is  provided  by  Plaut  et  al.  (1994).  The  method  is 
widely  used  to  identify  coherent  space-time  patterns  from  historical  observations  at 
regular  and  irregular  locations.  Unlike  traditional  spectrum  analysis,  where  the  basis 
functions  are  sine  and  cosines,  SSA  is  determined  from  estimates  of  the  lagged  cross¬ 
covariance  (space  and/or  time)  where'basis  functions  are  data-adaptive,  empirical  and 
orthogonal.  The  method  can  be  shown  to  be  optimal  in  the  sense  of  capturing  the 
maximum  variance  with  the  fewest  independent  components.  The  goal  here  is  to 
estimate  periodic  or  nearly  periodic  variance  components  in  P-T-R  and  relate  these  to 
hydrologic  conditions  across  the  Wasatch  Front.  A  brief  description  of  the  method 
follows. 


Theory 

We  first  define  an  L-dimensional  data  vector  Xu,  1  <1<L,  l<i<N,  where  l  is  channel 
number  representing  spatial  variables,  i  represents  the  time  index,  N  is  record  length 
and  the  mean  of  each  time  series  has  been  removed.  The  expansion  of  Xu  is  given  by: 

*u+,=  1^,1  (27) 

k-\ 

where  M is  the  window  length  or  embedding  dimension.  Thus,  X,i+I>  Xli+2, . ..,  X,i+M 

represent  shifted  time  series  at  channel  /  from  the  multichannel  record  Xu.  The 
coefficient  4*  is  called  the  space-time  principal  component  (ST-PC).  The  vector  £*, 
is  defined  as  a  space-time  empirical  orthogonal  function  (ST-EOF),  and  is  the  kth 
eigenvector  of  the  block-Toeplitz  matrix  Tx: 
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Tlt  is  the  MxM  lag  covariance  matrix  between  channel  l  and  channel  l  with  the 

element  The  least-biased  estimate  of  the  element  (7),.)^  in  row  j  and  column  f 
is  [Vautard  et  al.,1992]: 

Mr*  (29) 

The  ST-PCs  of  Tx  are  defined  as  the  orthogonal  projections  of  the  original  series  on 
the  ST-EOFs: 

A  =  XE.  (30) 

The  dimension  of  A  is  ( N-M+l)x(L+M ).  The  ST-PCs  can  be  thought  of  as  weighted 

moving  averages  of  the  original  time  series  X. 

We  are  mainly  interested  in  the  ST-PC’s  which  contribute  a  large  fraction  of  the 
variance  and  those  which  represent  oscillatory  components  of  the  signal  and  can  in 
some  way  be  distinguished  from  the  underlying  noise  in  the  record.  We  will  revisit 
the  question  of  resolving  significant  periodic  or  almost  periodic  components  and 
noise  later.  The  standard  approach  is  to  sort  the  eigenvalues  At  in  descending  order. 

The  plot  of  eigenvalue  A*  versus  k  represents  a  ranking  of  the  relative  variance 

contribution  for  each  ST-PC.  If  adjacent  eigenvalues  are  nearly  equal  and  their 
respective  eigenvectors  are  in  quadrature,  then  they  are  assumed  to  be  an  oscillatory 
pair  and  the  ST-EOF  pairs  are  either  periodic  or  almost  periodic.  The  oscillatory  pair 
would  represent  a  sine  and  cosine-pair  with  a  phase  shift.  We  apply  this  approach 
here,  inspecting  the  £*’ s  of  the  largest  eigenvalues  A*  from  which  the  period 

(frequency)  is  estimated  for  a  particular  channel  or  altitude. 

The  cutoff  for  “significant”  components  is  arbitrarily  determined  to  be  the 
value  of  k  where  a  flattening  or  change  in  slope  of  the  ranked  eigenvalues  is  observed 
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and  the  variance  contribution  is  small.  This  break  in  slope  is  sometimes  referred  to  as 
the  “noise  floor”.  Of  course  this  definition  is  problematic  since  ranking  of 
eigenvalues  by  size  has  no  ordering  with  respect  to  the  period  or  frequency  content  of 
the  data,  as  would  be  the  case  for  Fourier  spectral  estimates.  It  is  possible  that 
significant  oscillatory  pairs  may  not  be  detected,  especially  in  cases  where  the 
underlying  noise  spectrum  is  from  a  correlated  stochastic  process,  which  is  often  the 
case  in  hydrologic  data.  Next  we  examine  the  space-time  characteristics  of 
precipitation  -temperature  and  runoff  (P-T-R)  for  the  Wasatch  Front,  northern  Utah, 
using  multi-channel  SSA  and  the  qualitative  approach  to  detecting  oscillations 
described  above. 


State-Space  Dimension 

Multichannel  SSA  was  applied  separately  to  the  precipitation,  temperature  and 
runoff  (P-T-R)  data  sets.  The  concurrent  record  was  from  January  1948  to  December 
1990,  representing  43  years  of  data.  As  in  Fourier  spectrum  and  covariance  analysis, 
it  is  necessary  to  test  the  effect  of  the  time  embedding  (M)  or  lag  on  SSA  estimates. 
The  range  140^Vfel72  was  examined  and  M=160  seemed  to  provide  the  best 
resolution  across  the  range  of  oscillations  of  interest  (annual  to  decadal).  The  ranked 
eigenvalues  for  each  variable  is  shown  in  Figure  7.  From  the  graph,  the  near  equality 
of  some  pairs  of  leading  eigenvalues  is  noted.  These  pairs  often  (but  not  always) 
correspond  to  eigenvectors  that  are  periodic  or  nearly  periodic  and  in  quadrature.  The 
analysis  focuses  on  the  first  8-10  components,  up  to  the  point  where  the  slope  of  Xk 

versus  k  flattens  and  the  variance  contribution  is  small.  A  summary  of  oscillatory 
pairs  and  the  estimated  period  for  precipitation,  temperature  and  runoff  is  given  in 
Table  2. 

The  annual  oscillation  represents  the  leading  pair  of  eigenvalues  for  P-T-R  as  shown 
in  Table  2.  For  precipitation  the  annual  oscillation  represents  31%  of  the  total  space- 
time  variance.  Annual  harmonics  at  periods  0.5, 0.33  and  0.25  years  individually 
contribute  from  1.5- 1.9%  of  the  total  variance  for  P.  The  period  of  1 1  years  accounts 
for  1.7%  of  the  variance,  while  the  interannual  time  scale,  period  3  contributes  1.1%. 
In  the  case  of  temperature,  the  annual  oscillatory  pair  represents  93%  of  the  total 
variance  and  the  semi-annual  pair  contributes  an  additional  1%.  Interannual  and 
decadal  oscillations  are  all  quite  small,  with  the  interannual  component  contributing 
-0.1%. 

For  runoff  the  annual  oscillation  contributes  27%  of  the  total  variance  of  the  9 
stations,  while  the  second  largest  oscillatory  pair  is  at  the  1 1-year  period  with  10.5%. 
The  seasonal  harmonics  contribute  another  12%,  while  the  interannual  components 
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represent  2.6%  of  the  total  variance.  Overall,  we  can  say  that  the  variance 
contribution  for  P-T-R  at  all  altitudes  is  dominated  by  the  annual  oscillation  and 
harmonics,  precipitation  shows  significant  but  weak  interannual  and  decadal 
oscillations,  while  runoff  exhibits  strong  interannual  and  decadal  oscillations. 


Qrderfe 


Figure  7.  The  eigenvalue  spectrum  for  precipitation,  temperature  and  runoff  across  the  Wasatch 
Front,  northern  Utah.  The  "noise  floor"  gives  a  measure  of  the  statistical  dimension  of  the  hydrologic 
system. 
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Reconstruction  of  the  space-time  structure  of  P-T-R 


In  order  to  get  a  better  sense  of  how  well  the  first  few  eigenmodes  describe  the  P-T-R 
signal  above  the  noise  floor,  we  reconstruct  the  P-T-R  with  the  first  eight 
eigenvectors  and  examine  the  noise-free  response.  The  reconstructed  component  R 
associated  with  the  k[h  ST-EOF  at  time  i,  and  for  channel  /  is  given  by  [Plaut  et 
al.,1994]: 

i  m 

/£  =  -'LKjKj  when  1  <i  <M-1, 

1  M 

1  A/ 

Rik.  = -  Y  A*  £*  when  N-M  <i  <N  (31) 

A/  — j  +  1  ^  J  1,1 

IV  -  iTl  j=i-N+M 

where  the  range  of  k  is  from  1  to  LxM.  A  time  series  for  any  channel  can  be 

reconstructed  by  summing  over  the  dominant  components.  Figure  8  shows  the 
reconstructed  time  series  of  P-T-R  for  two  channels,  one  at  high  elevation  and  another 
at  a  low  elevation.  The  reconstructed  temperature  time  series  represents  95%  of  the 
total  variance  (4  ST-PC’s)  in  the  original  data.  For  precipitation,  seasonal  harmonics 
and  low-frequency  components  (10  ST-PC's)  explain  37%  of  the  total  variance,  while 
59%  of  the  total  variance  is  explained  for  the  runoff  data  (10  ST-PC’s). 


1  M 


when  M  <  i  <N-M+1 
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Brecipitslion: 


The  dynamic  behavior  of  the  hydrologic  system  can  be  examined  by  constructing  the 
time-  trajectory  or  parametric  curve  formed  by  using  the  reconstructed  state  variables 
as  coordinates.  Figure  9  illustrates  the  P-T-R  trajectory  at  three  elevations  of  Wasatch 
Range  for  the  reconstructed  time  series.  The  variance  components  were  the  same  as  in 
Figure  8  and  the  period  of  reconstruction  is  43  years.  At  the  highest  altitude  the  three 
dimensional  trajectory  exhibits  tightly  spaced  loops,  where  each  loop  represents  the 
annual  oscillation.  By  comparison,  at  low  elevation  the  phase  plane  is  smeared  out 
reflecting  the  introduction  of  the  low-frequency  components  and  deeper  groundwater 
discharge.  The  higher  elevations  are  clearly  dominated  by  seasonal  fluctuations  while 
the  low  elevations  contribute  the  longer  time  scales  imparted  by  the  characteristic 
hydrogeologic  conditions  controlling  runoff. 
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Figure  9.  A  phase-plane  plot  of  reconstructed  precipitation-temperature  -runoff  for  3 
elvations.  Note  the  low-frequency,  quasi-periodic  character  of  the  lowest  elevation 
where  the  long  time  scales  of  deep  groundwater  baseflow  to  streams  affects  the 
runoff. 


Figure  10.  The  hydrogeologic  conceptual  model  for  the  Wasatch  Front  showing  the  relation 
between  surface  and  subsurface  hydrology  and  demonstrating  the  importance  of  the  alluvial  fans  to 
supporting  low-frequency  contribution  to  streamflow  at  the  lowest  elevation.. 


Model  and  Parameter  Identification 

LDMS  applies  the  Genetic  Algorithm  (Holland,  )  to  the  problem  of  model  and  parameter 
identification.  There  are  many  system  identification  techniques  that  have  been  used  to  build  accurate 
mathematical  models  of  dynamic  systems.  These  are  generally  based  least  mean  squares  or 
maximum  likelihood  fitting  of  model  to  data,  and  most  are  some  sort  of  gradient-guided  local  search 
techniques  [Tan,  et  al.,  1995].  These  conventional  techniques  require  a  smooth  search  space  and 
easily  fail  in  obtaining  the  global  optimum  when  the  search  space  is  not  differentiable.  Furthermore, 
they  are  not  easily  be  applied  to  nonlinear  systems  and  are  not  valid  if  the  nonlinearities  of  the 
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system  are  not  differentiable  such  as  with  threshold  nonlinearities.  Genetic  Algorithms  (GAs)  are 
based  on  the  metaphor  of  natural  selection  and  Darwinian  evolution.  Genetic  Algorithms 
simultaneously  search  the  solutions  through  multiple  points  in  the  solution  space  and  often  are  found 
to  be  the  most  effective  method  of  finding  the  global  optimum. 

The  optimum  parameters  are  determined  in  such  a  way  that  the  sum  of  squares  of  the  difference 
between  the  actually  observed  and  the  modeled  values,  multiplied  by  a  number  that  measures  the 
degree  of  precision,  is  a  minimum.  This  function  is  called  the  fitness  function  is  given  by: 


In  our  case  q0  is  observed  runoff,  qm  is  modeled  runoff,  N  is  the  total  time  steps  of  the  simulation,  /  is 
the  time  index  and  a  is  a  positive  scaling  factor.  By  including  the  exponential  function,  the 
minimization  problem  changed  to  maximum  and  the  fitness  value  can  be  constrained  to  the  range  of 
(0,1]. 

In  order  to  examine  the  accuracy  of  each  models,  we  apply  several  forecast  accuracy  measures 
[Maidment,1993]  for  estimating  the  forecast  error.  The  measures  include  (1)  mean  absolute  error 
(MAE),  (2)  Relative  mean  absolute  error  (RMAE),  and  (3)  mean  square  error  (MSE) 
[Baidment,1993],  which  are  defined  by 

N  i=l 

RMAE  =  „ -  MAE  (33) 

Z«.(0 

1=1 

MSE  =  -|-  2[qji) -«„(()] 

/V  ,=1 

Field  Application:  West  Branch,  Susquehanna  River  Basin 

Modeling  the  rainfall-runoff  process  for  a  hydrologic  system  (river  basin,  drainage  basin,  catchment 
or  watershed)  has  been  a  central  problem  for  contemporary  hydrology  research.  The  problem  is 
made  difficult  because  of  the  heterogeneous  space-time  nature  of  hydrologic  systems.  Increasingly, 
complex  models  for  lareg-scale  hydrological  response  have  been  developed  using  tools  of  high 
performance  computation.  However,  as  has  been  pointed  out  by  Jakeman  and  Homberger  [1993],  it 
is  not  necessarily  true  that  spatially  detailed  complex  models  improve  the  model  results  and  in  fact 
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may  lead  to  over-parameterization.  That  is,  parameters  themselves  introduce  noise  or  model  errors 
and  limit  the  predictability  of  the  phenomenon. 

Geohvdrologic  Framework 

The  selected  research  site  is  the  Upper  West  Branch,  which  is  a  sub-basin  of  the  Susquehanna  River 
Basin.  Figure  2  describes  its  location.  This  river  basin  is  located  in  north  central  Pennsylvania  and 
has  an  area  of  14,710  km2.  Most  of  the  Upper  West  Branch  is  within  the  Appalachian  Plateau 
physiographic  province.  Only  a  very  small  portion  of  the  south  part  belongs  to  the  Appalachian 
Mountain  section  of  the  Valley  and  Ridge  physiographic  province.  The  regional  dip  of  the  bedrock 
rarely  exceeds  25  feet  per  mile.  Elongate,  gentle  folds  form  alternating  anticlines  and  synclines. 
These  gentle  warps  are  surface  expressions  of  small  displacements  along  deep-seated  faults. 
Deformation  of  strata  in  the  Appalachian  Plateaus  is  extremely  subde.  Glacial  erosion  and  deposition 
left  a  nearly  continuous  layer  of  till  over  upland  bedrock.  Colluvium  and  alluvial  deposits  mostly 
distribute  in  main  valleys.  In  this  region,  unconsolidated  glacial,  outwash,  colluvium  and  alluvial 
deposits  form  the  upper-layer  aquifer.  With  the  enough  thickness,  this  aquifer  is  mainly  productive. 
The  consolidated  sedimentary  rocks  mainly  are  shales,  sandstones,  siltstones  and  coal  seams,  which 
are  almost  flat-lying  to  gently  folded.  The  main  water-yielding  rocks  are  sandstones.  The  fractures, 
joints  and  pores  in  these  consolidated  rocks  form  the  bedrock  aquifers.  This  aquifer  is  constructed 
the  lower-layer  aquifer  in  the  Upper  West  Branch,  which  is  mostly  overlain  by  the  upper 
unconsolidated  aquifer. 

Figure  1 1.  depicts  the  hydrogeologic  conceptual  model  of  the  Upper  West  Branch.  The  water 
is  occurring  in  the  unconsolidated  aquifers  and  bedrock  aquifers.  Most  of  the  precipitation  directly 
recharges  into  the  upper  unconsolidated  layer.  In  the  upper  land,  since  the  till  is  thin  and  very 
permeable,  water  in  till  rapidly  discharges  to  colluvium  or  alluvial  aquifers.  In  lower  lands,  the 
precipitation  mainly  infiltrates  into  colluvium  and  alluvial  aquifers.  Then  groundwater  in  this  layer 
may  runoff  to  the  local  tributaries  or  streams.  This  forms  the  local  flow  system.  However,  a  small 
part  of  water  in  unconsolidated  aquifers  may  infiltrate  down  to  the  fractures  in  lower  bedrock 
aquifers.  Bedrock  aquifers  in  the  unglaciated  part  of  the  Upper  West  Branch  accept  less  direct 
recharge  from  the  precipitation  because  this  part  is  highly  dissected  and  much  of  the  area  is  near 
highly  sloping.  The  circulation  of  water  in  bedrock  aquifers  does  not  extends  no  more  than  a  few 
hundred  feet  below  the  land  surface  [Harlow  et  al,  1993]  while  the  permeability  of  bedrock  aquifers 
decreases  with  depth.  The  small  portion  of  the  water  in  unconsolidated  aquifer  may  moves  vertically 
downward  through  fractures,  then  horizontally  flows  through  the  consolidated  aquifers  and  finally 
discharges  to  the  streams.  Because  of  the  limited  flow  depth,  the  water  flow  paths  are  typically  not 
too  long,  commonly  extending  no  more  than  several  kilometers  in  their  longest  dimension  [Seaber, 
et  al.,  1988].  The  water  circulation  in  bedrock  aquifers  forms  the  intermediate  flow  system.  No 
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regional  flow  occurs  in  this  basin  because  the  thickness  of  the  regional  permeable  bedrock  aquifers 
is  less  than  several  hundred  feet. 
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Figure  11.  The  hydrogeologic  conceptual  model  for  a  typical  stream- reach  within  the  Upper  West 
Branch  of  the  Susquehanna  River.  Shallow  till,  colluvium  and  alluvium  form  the  near-surface 
groundwater  reservoir  while  the  deeper  subsurface  reservoir  is  dominanted  by  fracturing  in  various 
sedimentary  units. 


The  Model 

For  the  hydrogeologic  characteristics  of  the  Upper  West  Branch,  a  series  of  coupled  nonlinear 
integral-balance  model  equations  were  proposed  (Table  1).  This  serves  to  constrain  the  model  to  up 
to  two  state  variables  and  a  minimal  number  of  parameters.  The  six  models  of  the  Upper  West 
Branch  were  tested  and  parameters  estimated  with  the  Genetic  Algorithm  method  for  observed  daily 
precipitation,  temperature  and  streamflow  over  the  period  July  1985  to  June  1987.  The  modeled  and 
observed  data  are  depicted  in  Figure  12.  Comparison  of  the  results  of  the  linear  and  nonlinear 
models,  we  see  that  the  nonlinear  models  improve  the  peakflow  better  than  the  linear  models.  It  was 
also  determined  that  the  two-state  coupled  nonlinear  model  best  matches  the  observed  runoff.  The 
developed  models  also  were  verified  by  applying  the  calibrated  parameters  to  predict  the  runoff  in 
the  period  of  July  1989  to  June  1991  with  the  given  precipitation  and  temperature.  The  predicted 

30 


results  are  shown  in  Figure  13.  As  for  the  calibration  period,  the  model  results  for  the  two-state 
nonlinear  integral-balance  model  matches  the  observed  values  best. 


Model  Type 

Linear 

Nonlinear 

One-state 

s—  =  R-ay, 
dt  (16) 

qg  =  ccy; 

s^f  =  R-(ay  +  py2y 
dt  (17) 

qg=(ay+Py2); 

Two-state, 

uncoupled 

s^=(l-X)R-a,y2; 

dt 

q8  =  «i.vi  +  ovyo; 

(18) 

Sx^j-  =  M-(alyl+l3ly2ly, 
s2^=(l-m-(a2y2+p2y2y, 

at 

qs  =  (aji  +  ptf)  +(a2y2  +  j32y2 ); 
(19) 

Two-state,  coupled 

dy 

5i =  R-  ctxyx  - y(y,  -y0); 

dy,  ,  . 

j^-rCy, -*)-«*; 

qg  =  alyi  +  a1y2- 

(20) 

si—j-=R-(aiyi+  Pm)  -  r  0 i  -  %); 

J2  %  -3b)“  (^2  +/32>,2); 

at 

qs  =  (a,)i  +  Ay,2)  +(a2y2 + Piy\y 
(21) 

Table  2.  List  of  the  integral-balance  models  for  the  different  conceptualized  models  and  constitutive 
relationships  between  state  variables  and  output  flux  components. 
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Figure  12.  The  modeled  and  observed  runoff  for  the  Upper  West  Branch  for  each  of  the  6  models 
proposed  in  Table  ?. 
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Figure  13.  Forecast  for  the  Upper  West  Branch,  Susquehanna  River  using  all  6  models  and  the 
observed  precipitation  and  temperature  record. 


Accuracy  Analysis 


In  order  to  evaluate  the  accuracy  of  each  model,  we  used  several  forecast  accuracy  measures 
[Maidment,1993],  to  estimate  forecast  error.  The  measures  include  (1)  mean  absolute  error  (MAE), 
(2)  Relative  mean  absolute  error  (RMAE),  and  (3)  mean  square  error  (MSE)  [Baidment,1993], 
which  are  defined  by 

Ni= i 


RMAE  =■ 


MAE 


r=  1 

M5£  =  ^Z[?»(/)-9,(0]2 

Ar  i=i 
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The  lower  these  measures  are,  the  higher  the  prediction  accuracy  is.  Table  3  lists  the  estimations  of 
the  forecast  accuracy  analysis  for  the  identified  models. 


Model 

MAE 

RMAE 

MSE 

One-state  linear  system 

1.214 

0.3904 

3.9356 

Two-state  uncoupled  linear 
system 

1.2917 

0.4154 

4.3229 

Two-state  coupled 
linear  system 

1.1578 

0.3724 

3.541 

One-state  nonlinear  system 

1.046 

0.3364 

3.3625 

Two-state  uncoupled 
nonlinear  system 

1.0295 

0.3311 

3.1952 

Two-state  coupled 
nonlinear  system  . 

0.9764 

0.3140 

2.7560 

Table  3.  Forecast  accuracy  for  the  six  developed  models  of  the  Upper  West  Branch. 

We  note  that  the  two-state  coupled  nonlinear  version  of  the  model  represents  the  best  model 
performance  of  the  Upper  West  Branch. 

Figure  12  and  13  also  show  that  during  periods  of  drought  as  in  July,  August,  September  and 
October  in  both  1985  and  1989,  the  model  over-estimates  the  runoff.  We  propose  that  a  major  reason 
for  this  overestimation  is  that  unsaturated  groundwater  flow  is  not  included  in  the  initial  models.  In 
the  periods  of  drought,  the  soil  moisture  reservoir  may  totally  dominate  the  runoff  by  reducing 
recharge  to  saturated  groundwater  and  absorbing  all  precipitation  events.  When  rainfall  occurs,  the 
soil  moisture  first  has  to  meet  the  deficit  and  satisfy  soil  field  capacity  before  infiltration  and 
recharge  can  occur.  Future  improvements  to  the  model  will  explore  this  process. 
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